Ejercicios Modelos lineales

Ejercicio 1: Modelo regresión

Tenemos un dataset con riqueza de especies en las islas Galápagos. Nos preguntamos si hay un efecto de las variables: área de la isla, elevación máxima de la isla y distancia a la isla más próxima en la riqueza de especies.

Variables:

  • Area: área de la isla
  • Elevation: elevación máxima de la isla
  • Nearest: distancia isla más próxima
  • Species: riqueza de especies

Pasos que seguir

  1. Instalar y cargar el paquete faraway (data(gala))

  2. Representar gráficamente relación variable respuesta con las variables explicativas y correlación variables

  3. Ajustar modelo lineal

  4. Interpretación resultados modelo

  5. Comprobación supuestos del modelo

Instalar paquete y cargar los datos

library(faraway)
data(gala)
head(gala)
             Species Endemics  Area Elevation Nearest Scruz Adjacent
Baltra            58       23 25.09       346     0.6   0.6     1.84
Bartolome         31       21  1.24       109     0.6  26.3   572.33
Caldwell           3        3  0.21       114     2.8  58.7     0.78
Champion          25        9  0.10        46     1.9  47.4     0.18
Coamano            2        1  0.05        77     1.9   1.9   903.82
Daphne.Major      18       11  0.34       119     8.0   8.0     1.84

1. Análisis exploratorios

library(ggplot2)
library(patchwork)
plot1 <- ggplot(data = gala, aes(x = Area , y = Species)) +  
  geom_point(color = "blue", size = 1.5) +      
  labs(x = "Área de la isla",         
       y = "Riqueza de especies",                   
      title = "Relación entre riqueza de especies y área de la isla"   
  ) +
  theme_minimal()   
plot2 <- ggplot(data = gala, aes(x = Elevation , y = Species)) +  
  geom_point(color = "red", size = 1.5) +      
  labs(x = "Elevación de la isla",         
       y = "",                   
      title = "Relación entre riqueza de especies y elevación de la isla"   
  ) +
  theme_minimal()
plot3 <- ggplot(data = gala, aes(x = Nearest , y = Species)) +  
  geom_point(color = "darkgreen", size = 1.5) +      
  labs(x = "Distancia isla más próxima",         
       y = "",                   
      title = "Relación entre riqueza de especies y distancia isla más próxima"   
  ) +
  theme_minimal() 
plot1 + plot2 + plot3

2. Correlación variables explicativas

library(GGally)
ggcorr(gala[,c(1,3:5)], label = TRUE)

Parece relación entre área y elevación. Convertir a log podría ayudar.

Relación entre área y elevación. ¿Qué se podría hacer? Quitar una.

3. Ajuste modelo lineal e interpretación modelo

lm.gal0<-lm(Species~Area+Elevation+Nearest, gala)
anova(lm.gal0)
Analysis of Variance Table

Response: Species
          Df Sum Sq Mean Sq F value    Pr(>F)    
Area       1 145470  145470 22.2592 7.076e-05 ***
Elevation  1  65664   65664 10.0476  0.003882 ** 
Nearest    1     29      29  0.0045  0.947179    
Residuals 26 169918    6535                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. Comprobación supuestos del modelo

par(mfcol=c(2,2))
plot(lm.gal0)

Ausencia de homocedasticidad, normalidad y linealidad. Y valores atípicos. No sirve modelo.

Factor de inflación de la varianza (VIF)

Cómo incrementa la varianza de un coeficiente estimado del modelo de regresión cuando los predictores están correlacionados entre sí

library(car)
vif(lm.gal0)
     Area Elevation   Nearest 
 2.373470  2.344460  1.025189

Transformación variable respuesta

gala$log_Sp <- log(gala$Species)
plot1 <- ggplot(data = gala, aes(x = Area , y = log_Sp)) +  
  geom_point(color = "blue", size = 1.5) +      
  labs(x = "Área de la isla",         
       y = "Riqueza de especies",                   
      title = "Relación entre riqueza de especies y área de la isla"   
  ) +
  theme_minimal()   
plot2 <- ggplot(data = gala, aes(x = Elevation , y = log_Sp)) +  
  geom_point(color = "red", size = 1.5) +      
  labs(x = "Elevación de la isla",         
       y = "",                   
      title = "Relación entre riqueza de especies y elevación de la isla"   
  ) +
  theme_minimal()
plot3 <- ggplot(data = gala, aes(x = Nearest , y = log_Sp)) +  
  geom_point(color = "darkgreen", size = 1.5) +      
  labs(x = "Distancia isla más próxima",         
       y = "",                   
      title = "Relación entre riqueza de especies y distancia isla más próxima"   
  ) +
  theme_minimal() 
plot1 + plot2 + plot3

Sigue sin ser lineal la relación.

Transformación variables predictoras

gala$log_area <- log(gala$Area)
gala$log_elevacion <- log(gala$Elevation)
gala$log_near <- log(gala$Nearest)

plot1 <- ggplot(data = gala, aes(x = log_area , y = log_Sp)) +  
  geom_point(color = "blue", size = 2) +      
  labs(x = "Área de la isla",         
       y = "Riqueza de especies",                   
      title = "Relación entre riqueza de especies y área de la isla"   
  ) +
  theme_minimal()   
plot2 <- ggplot(data = gala, aes(x = log_elevacion , y = log_Sp)) +  
  geom_point(color = "red", size = 2) +      
  labs(x = "Elevación de la isla",         
       y = "",                   
      title = "Relación entre riqueza de especies y elevación de la isla"   
  ) +
  theme_minimal()
plot3 <- ggplot(data = gala, aes(x = log_near , y = log_Sp)) +  
  geom_point(color = "darkgreen", size = 2) +      
  labs(x = "Distancia isla más próxima",         
       y = "",                   
      title = "Relación entre riqueza de especies y distancia isla más próxima"   
  ) +
  theme_minimal() 
plot1 + plot2 + plot3

Ajuste modelo y comprobación supuestos

lm.gal1 <- lm(log_Sp~log_area+log_elevacion+log_near, gala)
par(mfcol=c(2,2))
plot(lm.gal1)

¿Jerarquización variables?

Anova(lm.gal1,type="III")
Anova Table (Type III tests)

Response: log_Sp
               Sum Sq Df F value    Pr(>F)    
(Intercept)    5.4127  1  9.0225  0.005832 ** 
log_area      15.2592  1 25.4356 2.997e-05 ***
log_elevacion  0.6725  1  1.1210  0.299433    
log_near       1.0009  1  1.6684  0.207841    
Residuals     15.5978 26                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Problemas colinealidad

vif(lm.gal1)
     log_area log_elevacion      log_near 
     5.480281      5.461635      1.008644

Valor muy grande. Decidir qué variable eliminar

lm.gal2 <- lm(log_Sp~log_area+log_near,gala)
Anova(lm.gal2,type="III")
Anova Table (Type III tests)

Response: log_Sp
             Sum Sq Df  F value    Pr(>F)    
(Intercept) 156.581  1 259.8405 2.225e-15 ***
log_area     54.502  1  90.4448 4.123e-10 ***
log_near      0.948  1   1.5731    0.2205    
Residuals    16.270 27                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Quitar distancia

lm.gal3 <- lm(log_Sp~log_area,gala)
Anova(lm.gal3,type="III")
Anova Table (Type III tests)

Response: log_Sp
             Sum Sq Df F value    Pr(>F)    
(Intercept) 210.109  1 341.676 < 2.2e-16 ***
log_area     53.668  1  87.274  4.23e-10 ***
Residuals    17.218 28                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm.gal3)

Call:
lm(formula = log_Sp ~ log_area, data = gala)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5442 -0.4001  0.0941  0.5449  1.3752 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.9037     0.1571  18.484  < 2e-16 ***
log_area      0.3886     0.0416   9.342 4.23e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7842 on 28 degrees of freedom
Multiple R-squared:  0.7571,    Adjusted R-squared:  0.7484 
F-statistic: 87.27 on 1 and 28 DF,  p-value: 4.23e-10

log(Species) = 2.9 + 0.39 log(Area)

Gráfico de predicciones

library(broom)
library(tidyverse)
library(dplyr)
preds <- augment(lm.gal3, type.predict = "response", se_fit = TRUE)
preds <- preds %>%
  mutate(lower = .fitted - 1.96 * .se.fit,upper = .fitted + 1.96 * .se.fit)
  
ggplot(data = preds, aes(x = log_area, y = log_Sp)) +  
  geom_point(color = "blue", size = 2) +     
  geom_line(aes(y = .fitted)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
  labs(x = "Área de la isla",         
       y = "Riqueza de especies",                   
      title = "Predicciones modelo de regresión"   
  ) +
  theme_minimal()

Ejercicio 2: Modelo factorial

Tenemos un dataset con biomasa leñosa en un área mediterránea. Nos preguntamos si hay un efecto de la variable cobertura arbórea considerando dos niveles de sequía estival.

Variables:

  • biomasa: biomasa leñosa
  • cobertura: cobertura arbórea
  • sequia: sequía o lluvia estival

Pasos que seguir

  1. Cargar los datos (cobertura.csv)

  2. Representar gráficamente relación variable respuesta con las variables explicativas

  3. Ajustar modelo lineal

  4. Interpretación resultados modelo

  5. Comprobación supuestos del modelo

Instalar paquete y cargar los datos

library(here)
cobertura <- read.csv(here("data/cobertura.csv"))
head(cobertura)
  X  biomasa cobertura sequia
1 1 44.99793  18.55938   seco
2 2 24.05575  24.14701   seco
3 3 62.48798  68.89674   seco
4 4 47.19165  30.01684   seco
5 5 52.97670  31.25460   seco
6 6 54.23581  38.79990   seco

1. Análisis exploratorios

library(ggplot2)
library(patchwork)
ggplot(cobertura, aes(x = biomasa, y = cobertura)) +
  geom_point(color = "steelblue", size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  facet_wrap(~ sequia) +
  labs(
    title = "Relación biomasa leñosa y cobertura arbórea",
    x = "Cobertura arbórea",
    y = "Biomasa leñosa"
  ) +
  theme_minimal()

Ajuste modelo de regresión

cobertura <- cobertura %>%
  mutate(sequia = as.factor(sequia))
lm_bio1 <- lm(biomasa ~ cobertura*sequia, data = cobertura)

Comprobación de hipótesis

anova(lm_bio1)
Analysis of Variance Table

Response: biomasa
                  Df Sum Sq Mean Sq  F value    Pr(>F)    
cobertura          1 169328  169328 1149.130 < 2.2e-16 ***
sequia             1  28187   28187  191.286 < 2.2e-16 ***
cobertura:sequia   1   7224    7224   49.023 3.936e-11 ***
Residuals        196  28881     147                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Supuestos del modelo

par(mfrow=c(2,2))
plot(lm_bio1)

Interpretación del modelo

summary(lm_bio1)

Call:
lm(formula = biomasa ~ cobertura * sequia, data = cobertura)

Residuals:
    Min      1Q  Median      3Q     Max 
-39.708  -6.943   0.206   7.332  51.094 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          25.91031    3.26354   7.939 1.54e-13 ***
cobertura             1.33964    0.06143  21.808  < 2e-16 ***
sequiaseco           -1.78935    4.20599  -0.425    0.671    
cobertura:sequiaseco -0.72258    0.10320  -7.002 3.94e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.14 on 196 degrees of freedom
Multiple R-squared:  0.8764,    Adjusted R-squared:  0.8745 
F-statistic: 463.1 on 3 and 196 DF,  p-value: < 2.2e-16

yi = 25.91 - 1.79 seco + (1.34 - 0.72 seco) cobertura

Gráfica de predicciones

library(broom)
newdata <- expand.grid(
  cobertura = seq(min(cobertura$cobertura), max(cobertura$cobertura), length.out = 100),
  sequia = levels(cobertura$sequia)
)
predicciones_grid <- augment(lm_bio1, newdata = newdata)

ggplot(cobertura, aes(x = cobertura, y = biomasa, color = sequia)) +
  geom_point(alpha = 0.5) +
  geom_line(data = predicciones_grid, aes(x= cobertura, y = .fitted, group = sequia, color = sequia), size = 1.2) +
  scale_color_manual(values = c("#B3E5FC", "#F8BBD0")) +
  labs(
    title = "ANCOVA: predicciones por nivel de sequía",
    x = "Cobertura",
    y = "Biomasa",
    color = "Sequía"
  ) +
  theme_minimal()